Surface emitter
In [1]:
Copied!
# Useful for debugging
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'
import os
import yaml
import numpy as np
# Useful for debugging
%load_ext autoreload
%autoreload 2
%config InlineBackend.figure_format = 'retina'
import os
import yaml
import numpy as np
In [2]:
Copied!
from distgen import Generator
from distgen import Generator
In [3]:
Copied!
input_yaml = """
n_particle: 100000
species: electron
start:
type: free
random:
type: hammersley
total_charge:
units: pC
value: 10
r_dist:
max_r:
units: mm
value: 1
type: radial_uniform
p_dist:
sigma_p:
value: 1
units: eV/c
avg_p:
value: 100
units: eV/c
type: gaussian
"""
input_yaml = """
n_particle: 100000
species: electron
start:
type: free
random:
type: hammersley
total_charge:
units: pC
value: 10
r_dist:
max_r:
units: mm
value: 1
type: radial_uniform
p_dist:
sigma_p:
value: 1
units: eV/c
avg_p:
value: 100
units: eV/c
type: gaussian
"""
In [4]:
Copied!
gen = Generator(input_yaml, verbose=True)
gen = Generator(input_yaml, verbose=True)
In [5]:
Copied!
P = gen.run()
P = gen.run()
Distribution format: None
Warning: no output file specified, defaulting to "None".
Output file: None
Creating beam distribution....
Beam starting from: free
Total charge: 10 pC.
Number of macroparticles: 100000.
Assuming cylindrical symmetry...
r distribution: radial uniform
min_r = 0 mm, max_r = 1 mm
theta distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
p distribution: Gaussian
avg_p = 100 eV/c, sigma_p = 1.000 eV/c
azimuthal angle distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
polar angle distribution: uniform phi
min_phi = 0 rad, max_phi = 3.14159 rad
Shifting avg_x = -4.85224E-08 mm -> 0 mm
Scaling sigma_x = 0.499991 mm -> 0.5 mm
Shifting avg_y = 1.44675E-06 mm -> 0 mm
Scaling sigma_y = 0.49999 mm -> 0.5 mm
Shifting avg_px = -3.3074E-05 eV/c -> -6.12323E-15 eV/c
Scaling sigma_px = 57.7382 eV/c -> 57.7379 eV/c
Shifting avg_py = 1.54491E-06 eV/c -> 0 eV/c
Scaling sigma_py = 57.7382 eV/c -> 57.7379 eV/c
Shifting avg_pz = -0.000118434 eV/c -> 6.12323E-15 eV/c
Scaling sigma_pz = 57.7372 eV/c -> 57.7379 eV/c
...done. Time Elapsed: 107.992 ms.
Created particles in .particles:
ParticleGroup with 100000 particles with total charge 1.0000000000000003e-11 C
In [6]:
Copied!
P.plot("p")
P.plot("p")
In [7]:
Copied!
P.plot("x", "px")
P.plot("x", "px")
In [8]:
Copied!
P.plot("pz")
P.plot("pz")
In [9]:
Copied!
P.plot("px")
P.plot("px")
In [10]:
Copied!
P.plot("py")
P.plot("py")
In [11]:
Copied!
from matplotlib import pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
px, py, pz = P["px"], P["py"], P["pz"]
p = np.sqrt(px**2 + py**2 + pz**2)
theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2
x = px / p
y = py / p
z = pz / p
ax.scatter(x[::200], y[::200], z[::200], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")
from matplotlib import pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
px, py, pz = P["px"], P["py"], P["pz"]
p = np.sqrt(px**2 + py**2 + pz**2)
theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2
x = px / p
y = py / p
z = pz / p
ax.scatter(x[::200], y[::200], z[::200], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")
Out[11]:
Text(0.5, 0, '$\\hat{p}_z$')
In [12]:
Copied!
theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2
plt.plot(pcs, hist);
theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2
plt.plot(pcs, hist);
In [13]:
Copied!
rns = np.random.random(100000)
rns = np.random.random(100000)
In [14]:
Copied!
pa = 0
pb = np.pi
Ca = np.cos(pa)
Cb = np.cos(pb)
pa = 0
pb = np.pi
Ca = np.cos(pa)
Cb = np.cos(pb)
In [15]:
Copied!
rho = 1 / (Ca - Cb)
rho = 1 / (Ca - Cb)
In [16]:
Copied!
phis = np.linspace(pa, pb, 1000)
phis = np.linspace(pa, pb, 1000)
In [17]:
Copied!
plt.plot(phis, rho * np.sin(phis));
plt.plot(phis, rho * np.sin(phis));
In [18]:
Copied!
np.trapezoid(rho * np.sin(phis), phis)
np.trapezoid(rho * np.sin(phis), phis)
Out[18]:
np.float64(0.999999175885426)
In [19]:
Copied!
cdf = (Ca - np.cos(phis)) * rho
cdf = (Ca - np.cos(phis)) * rho
In [20]:
Copied!
from scipy.integrate import cumulative_trapezoid as cumtrapz
from scipy.integrate import cumulative_trapezoid as cumtrapz
In [21]:
Copied!
plt.plot(phis, cdf, phis, cumtrapz(rho * np.sin(phis), phis, initial=0));
plt.plot(phis, cdf, phis, cumtrapz(rho * np.sin(phis), phis, initial=0));
In [22]:
Copied!
ps = np.arccos(Ca - rns * (Ca - Cb))
ps = np.arccos(Ca - rns * (Ca - Cb))
In [23]:
Copied!
hist, pedges = np.histogram(ps, bins=30, density=True)
hist, pedges = np.histogram(ps, bins=30, density=True)
In [24]:
Copied!
pcs = (pedges[1:] + pedges[:-1]) / 2
pcs = (pedges[1:] + pedges[:-1]) / 2
In [25]:
Copied!
plt.plot(pcs, hist, phis, rho * np.sin(phis));
plt.plot(pcs, hist, phis, rho * np.sin(phis));
In [26]:
Copied!
gen = Generator("data/beer.can.in.yaml", verbose=1)
# print(gen)
gen = Generator("data/beer.can.in.yaml", verbose=1)
# print(gen)
In [27]:
Copied!
P1 = gen.run()
P1 = gen.run()
Distribution format: gpt
Output file: beer.can.out.txt
Creating beam distribution....
Beam starting from: cathode
Total charge: 10 pC.
Number of macroparticles: 200000.
Assuming cylindrical symmetry...
r distribution: radial uniform
min_r = 0 mm, max_r = 2 mm
theta distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
t distribution: uniform
min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t: 1.1547 ps
px distribution: Gaussian
avg_px = 0 eV/c, sigma_px = 276.857 eV/c
py distribution: Gaussian
avg_py = 0 eV/c, sigma_py = 276.857 eV/c
pz distribution: Gaussian
avg_pz = 0 eV/c, sigma_pz = 276.857 eV/c
Shifting avg_x = -2.10539E-05 mm -> 0 mm
Scaling sigma_x = 0.999982 mm -> 1 mm
Shifting avg_y = -1.97648E-06 mm -> 0 mm
Scaling sigma_y = 0.999998 mm -> 1 mm
Shifting avg_px = -0.0212516 eV/c -> 0 eV/c
Scaling sigma_px = 276.849 eV/c -> 276.857 eV/c
Shifting avg_py = -0.0259627 eV/c -> 0 eV/c
Scaling sigma_py = 276.844 eV/c -> 276.857 eV/c
Shifting avg_pz = -0.0372953 eV/c -> 0 eV/c
Scaling sigma_pz = 276.843 eV/c -> 276.857 eV/c
Shifting avg_t = -4.90022E-05 ps -> 0 ps
Scaling sigma_t = 1.1547 ps -> 1.1547 ps
Cathode start: fixing pz momenta to forward hemisphere
avg_pz -> 220.904 eV/c, sigma_pz -> 166.887 eV/c
...done. Time Elapsed: 197.597 ms.
Created particles in .particles: ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C
In [28]:
Copied!
P1.plot("x", "px")
P1.plot("x", "px")
In [29]:
Copied!
P1.plot("p")
P1.plot("p")
In [30]:
Copied!
P1.plot("pz")
P1.plot("pz")
In [31]:
Copied!
P1["sigma_px"], P1["sigma_py"], P1["sigma_pz"]
P1["sigma_px"], P1["sigma_py"], P1["sigma_pz"]
Out[31]:
(np.float64(276.8570795554991), np.float64(276.8570795554991), np.float64(166.88683722575215))
In [32]:
Copied!
gen = Generator("data/maxwell_boltzmann.beer.can.in.yaml", verbose=1)
gen = Generator("data/maxwell_boltzmann.beer.can.in.yaml", verbose=1)
In [33]:
Copied!
P2 = gen.run()
P2 = gen.run()
Distribution format: gpt
Output file: beer.can.out.txt
Creating beam distribution....
Beam starting from: cathode
Total charge: 10 pC.
Number of macroparticles: 200000.
Assuming cylindrical symmetry...
r distribution: radial uniform
min_r = 0 mm, max_r = 2 mm
theta distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
p distribution: Maxwell-Boltzmann
p scale = 276.857 eV/c
azimuthal angle distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
polar angle distribution: uniform phi
min_phi = 0 rad, max_phi = 3.14159 rad
t distribution: uniform
min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t: 1.1547 ps
Shifting avg_x = 1.63849E-05 mm -> 0 mm
Scaling sigma_x = 0.999984 mm -> 1 mm
Shifting avg_y = -2.77139E-06 mm -> 0 mm
Scaling sigma_y = 0.999996 mm -> 1 mm
Shifting avg_px = -0.000954618 eV/c -> -2.70524E-14 eV/c
Scaling sigma_px = 276.853 eV/c -> 276.857 eV/rad⁰⋅⁵/c
Shifting avg_py = 0.000856608 eV/c -> 0 eV/c
Scaling sigma_py = 276.848 eV/c -> 276.857 eV/rad⁰⋅⁵/c
Shifting avg_pz = -0.0142766 eV/c -> 2.70524E-14 eV/c
Scaling sigma_pz = 276.833 eV/c -> 276.857 eV/rad⁰⋅⁵/c
Shifting avg_t = -4.90022E-05 ps -> 0 ps
Scaling sigma_t = 1.1547 ps -> 1.1547 ps
Cathode start: fixing pz momenta to forward hemisphere
avg_pz -> 220.909 eV/c, sigma_pz -> 166.881 eV/c
...done. Time Elapsed: 197.52 ms.
Created particles in .particles: ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C
In [34]:
Copied!
P2["sigma_px"], P2["sigma_py"], P2["sigma_pz"]
P2["sigma_px"], P2["sigma_py"], P2["sigma_pz"]
Out[34]:
(np.float64(276.85707955549907), np.float64(276.85707955549907), np.float64(166.88083659251544))
In [35]:
Copied!
gen = Generator("data/maxwell_boltzmann_KE.beer.can.in.yaml", verbose=1)
gen = Generator("data/maxwell_boltzmann_KE.beer.can.in.yaml", verbose=1)
In [36]:
Copied!
P3 = gen.run()
P3 = gen.run()
Distribution format: gpt Output file: beer.can.out.txt Creating beam distribution.... Beam starting from: cathode Total charge: 10 pC. Number of macroparticles: 200000. Assuming cylindrical symmetry...
r distribution: radial uniform
min_r = 0 mm, max_r = 2 mm
theta distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
KE distribution: Maxwell-Boltzmann Energy
kT = 150 meV
azimuthal angle distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
polar angle distribution: uniform phi
min_phi = 0 rad, max_phi = 3.14159 rad
t distribution: uniform
min_t = -2 ps, max_t = 2 ps, avg_t = 0 ps, sigma_t: 1.1547 ps
Shifting avg_x = 1.63849E-05 mm -> 0 mm
Scaling sigma_x = 0.999984 mm -> 1 mm
Shifting avg_y = -2.77139E-06 mm -> 0 mm
Scaling sigma_y = 0.999996 mm -> 1 mm
Shifting avg_px = -1.75487E-09 meV·s/m -> -9.02058E-20 meV·s/m
Scaling sigma_px = 0.000922998 meV·s/m -> 0.000922972 meV·s/m
Shifting avg_py = 1.21968E-09 meV·s/m -> 0 meV·s/m
Scaling sigma_py = 0.000922978 meV·s/m -> 0.000922972 meV·s/m
Shifting avg_pz = -4.48933E-08 meV·s/m -> 9.02058E-20 meV·s/m
Scaling sigma_pz = 0.00092294 meV·s/m -> 0.000922972 meV·s/m
Shifting avg_t = -4.90022E-05 ps -> 0 ps
Scaling sigma_t = 1.1547 ps -> 1.1547 ps
Cathode start: fixing pz momenta to forward hemisphere
avg_pz -> 0.0007366 meV·s/m, sigma_pz -> 0.000556145 meV·s/m
...done. Time Elapsed: 205.344 ms.
Created particles in .particles:
ParticleGroup with 200000 particles with total charge 1.0000000000000003e-11 C
In [37]:
Copied!
P3.plot("kinetic_energy")
P3.plot("kinetic_energy")
In [38]:
Copied!
P3.plot("x", "px")
P3.plot("x", "px")
In [39]:
Copied!
P3["sigma_px"], P3["sigma_py"], P3["sigma_pz"]
P3["sigma_px"], P3["sigma_py"], P3["sigma_pz"]
Out[39]:
(np.float64(276.70000872169317), np.float64(276.70000872169317), np.float64(166.7279863589753))
In [40]:
Copied!
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
px, py, pz = P3["px"], P3["py"], P3["pz"]
p = np.sqrt(px**2 + py**2 + pz**2)
theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2
x = px / p
y = py / p
z = pz / p
ax.scatter(x[::500], y[::500], z[::500], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
px, py, pz = P3["px"], P3["py"], P3["pz"]
p = np.sqrt(px**2 + py**2 + pz**2)
theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2
x = px / p
y = py / p
z = pz / p
ax.scatter(x[::500], y[::500], z[::500], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")
Out[40]:
Text(0.5, 0, '$\\hat{p}_z$')
In [41]:
Copied!
KE = np.linspace(0, 150, 10000)
E0 = 0
A1 = (0.8,)
m1 = 8
A2 = 0.1
m2 = 90
PKE = A1 * np.exp(-np.abs(KE - E0) / m1) + A2 * np.exp(-np.abs(KE - E0) / m2)
PKE = PKE / np.trapz(PKE, KE) # Numerically intergate to normalize
plt.plot(KE, PKE)
plt.xlabel("KE (eV)")
plt.ylabel("$\\rho(KE)$ (1/eV)");
KE = np.linspace(0, 150, 10000)
E0 = 0
A1 = (0.8,)
m1 = 8
A2 = 0.1
m2 = 90
PKE = A1 * np.exp(-np.abs(KE - E0) / m1) + A2 * np.exp(-np.abs(KE - E0) / m2)
PKE = PKE / np.trapz(PKE, KE) # Numerically intergate to normalize
plt.plot(KE, PKE)
plt.xlabel("KE (eV)")
plt.ylabel("$\\rho(KE)$ (1/eV)");
/tmp/ipykernel_2270/3289071326.py:12: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`. PKE = PKE / np.trapz(PKE, KE) # Numerically intergate to normalize
In [42]:
Copied!
dat = np.zeros((len(KE), 2))
dat[:, 0], dat[:, 1] = KE, PKE
np.savetxt("KEdist.txt", dat, header="KE PKE", comments="")
dat = np.zeros((len(KE), 2))
dat[:, 0], dat[:, 1] = KE, PKE
np.savetxt("KEdist.txt", dat, header="KE PKE", comments="")
In [43]:
Copied!
q = 1000 * 1.60217663e-19
print(q)
input_yaml = """
n_particle: 1000
species: electron
start:
type: cathode
random:
type: hammersley
total_charge:
units: C
value: 1.60217663e-17
r_dist:
max_r:
units: mm
value: 14.6
type: radial_uniform
KE_dist:
file: KEdist.txt
units: eV
type: file1d
transforms:
sx:
avg_x:
units: millimeter
value: -50
type: set_avg x
sy:
avg_y:
units: millimeter
value: 50
type: set_avg y
sz:
avg_z:
units: millimeter
value: 65
type: set_avg z
output:
file: None
"""
q = 1000 * 1.60217663e-19
print(q)
input_yaml = """
n_particle: 1000
species: electron
start:
type: cathode
random:
type: hammersley
total_charge:
units: C
value: 1.60217663e-17
r_dist:
max_r:
units: mm
value: 14.6
type: radial_uniform
KE_dist:
file: KEdist.txt
units: eV
type: file1d
transforms:
sx:
avg_x:
units: millimeter
value: -50
type: set_avg x
sy:
avg_y:
units: millimeter
value: 50
type: set_avg y
sz:
avg_z:
units: millimeter
value: 65
type: set_avg z
output:
file: None
"""
1.6021766299999998e-16
In [44]:
Copied!
inputs = yaml.safe_load(input_yaml)
n = 1000
inputs["n_particle"] = n
inputs["output"]["file"] = "test_ion_writer.ion"
inputs["output"]["type"] = "simion"
inputs["total_charge"]["value"] = n * 1.60217663e-19
inputs["transforms"]["sx"]["avg_x"]["value"] = -50
inputs["transforms"]["sy"]["avg_y"]["value"] = 50
inputs["transforms"]["sz"]["avg_z"]["value"] = 55
inputs = yaml.safe_load(input_yaml)
n = 1000
inputs["n_particle"] = n
inputs["output"]["file"] = "test_ion_writer.ion"
inputs["output"]["type"] = "simion"
inputs["total_charge"]["value"] = n * 1.60217663e-19
inputs["transforms"]["sx"]["avg_x"]["value"] = -50
inputs["transforms"]["sy"]["avg_y"]["value"] = 50
inputs["transforms"]["sz"]["avg_z"]["value"] = 55
In [45]:
Copied!
gen = Generator(inputs, verbose=1)
gen = Generator(inputs, verbose=1)
In [46]:
Copied!
P = gen.run()
P = gen.run()
Distribution format: simion
Output file: /home/runner/work/distgen/distgen/docs/examples/test_ion_writer.ion
Creating beam distribution....
Beam starting from: cathode
Total charge: 1.60218E-16 C.
Number of macroparticles: 1000.
Assuming cylindrical symmetry...
r distribution: radial uniform
min_r = 0 mm, max_r = 14.6 mm
theta distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
KE distribution: KE-distribution file: "/home/runner/work/distgen/distgen/docs/examples/KEdist.txt"
azimuthal angle distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
polar angle distribution: uniform phi
min_phi = 0 rad, max_phi = 3.14159 rad
Shifting avg_x = 0.00564031 mm -> 0 mm
Scaling sigma_x = 7.29551 mm -> 7.3 mm
Shifting avg_y = 0.00120927 mm -> 0 mm
Scaling sigma_y = 7.29741 mm -> 7.3 mm
Shifting avg_px = -6.09774E-09 eV·s/m -> -9.89015E-22 eV·s/m
Scaling sigma_px = 1.11822E-05 eV·s/m -> 1.11773E-05 eV·s/m
Shifting avg_py = -2.50205E-09 eV·s/m -> 0 eV·s/m
Scaling sigma_py = 1.11987E-05 eV·s/m -> 1.11773E-05 eV·s/m
Shifting avg_pz = -2.52855E-08 eV·s/m -> 9.89015E-22 eV·s/m
Scaling sigma_pz = 1.11509E-05 eV·s/m -> 1.11773E-05 eV·s/m
Cathode start: fixing pz momenta to forward hemisphere
avg_pz -> 8.07962E-06 eV·s/m, sigma_pz -> 7.72348E-06 eV·s/m
Applying user supplied transform: "sx" = set_avg x...
Setting avg_x -> -50 mm.
Applying user supplied transform: "sy" = set_avg y...
Setting avg_y -> 50 mm.
Applying user supplied transform: "sz" = set_avg z...
Setting avg_z -> 55 mm.
...done. Time Elapsed: 26.6883 ms.
Created particles in .particles:
ParticleGroup with 1000 particles with total charge 1.6021766300000003e-16 C
In [47]:
Copied!
P.plot("x", "y", figsize=(5, 5))
P.plot("x", "y", figsize=(5, 5))
In [48]:
Copied!
P.plot("x", "px", figsize=(5, 5))
P.plot("x", "px", figsize=(5, 5))
In [49]:
Copied!
P.plot("kinetic_energy")
P.plot("kinetic_energy")
In [50]:
Copied!
from distgen.writers import writer
writer("simion", gen.beam(), "test_ion_writer.ion", params={"color": 0})
from distgen.writers import writer
writer("simion", gen.beam(), "test_ion_writer.ion", params={"color": 0})
Distribution format: simion
Output file: /home/runner/work/distgen/distgen/docs/examples/test_ion_writer.ion
Creating beam distribution....
Beam starting from: cathode
Total charge: 1.60218E-16 C.
Number of macroparticles: 1000.
Assuming cylindrical symmetry...
r distribution: radial uniform
min_r = 0 mm, max_r = 14.6 mm
theta distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
KE distribution: KE-distribution file: "/home/runner/work/distgen/distgen/docs/examples/KEdist.txt"
azimuthal angle distribution: uniform theta
min_theta = 0 rad, max_theta = 6.28319 rad
polar angle distribution: uniform phi
min_phi = 0 rad, max_phi = 3.14159 rad
Shifting avg_x = 0.00564031 mm -> 0 mm
Scaling sigma_x = 7.29551 mm -> 7.3 mm
Shifting avg_y = 0.00120927 mm -> 0 mm
Scaling sigma_y = 7.29741 mm -> 7.3 mm
Shifting avg_px = -6.09774E-09 eV·s/m -> -9.89015E-22 eV·s/m
Scaling sigma_px = 1.11822E-05 eV·s/m -> 1.11773E-05 eV·s/m
Shifting avg_py = -2.50205E-09 eV·s/m -> 0 eV·s/m
Scaling sigma_py = 1.11987E-05 eV·s/m -> 1.11773E-05 eV·s/m
Shifting avg_pz = -2.52855E-08 eV·s/m -> 9.89015E-22 eV·s/m
Scaling sigma_pz = 1.11509E-05 eV·s/m -> 1.11773E-05 eV·s/m
Cathode start: fixing pz momenta to forward hemisphere
avg_pz -> 8.07962E-06 eV·s/m, sigma_pz -> 7.72348E-06 eV·s/m
Applying user supplied transform: "sx" = set_avg x...
Setting avg_x -> -50 mm.
Applying user supplied transform: "sy" = set_avg y...
Setting avg_y -> 50 mm.
Applying user supplied transform: "sz" = set_avg z...
Setting avg_z -> 55 mm.
...done. Time Elapsed: 28.3866 ms.
In [51]:
Copied!
from scipy.constants import physical_constants
mc2 = 1e6 * physical_constants["electron mass energy equivalent in MeV"][0]
e_ = physical_constants["elementary charge"][0]
me = physical_constants["electron mass in u"][0]
def particle_group_to_SIMION(P, filename, color=0):
header = ";0"
simion_params = [
"TOB",
"MASS",
"CHARGE",
"X",
"Y",
"Z",
"AZ",
"EL",
"KE",
"CWF",
"COLOR",
]
# simion_units = {
# "TOB": "usec",
# "MASS": "amu",
# "CHARGE": "e",
# "X": "mm",
# "Y": "mm",
# "Z": "mm",
# "AZ": "deg",
# "EL": "deg",
# "CWF": "",
# "COLOR": "",
# }
data = np.zeros((len(P), len(simion_params)))
data[:, simion_params.index("TOB")] = P.t * 1e6 # [P.t] = sec, convert to usec
if P.species == "electron":
data[:, simion_params.index("MASS")] = np.full(len(P), me)
data[:, simion_params.index("CHARGE")] = np.full(len(P), -1)
else:
raise ValueError(f"Species {P.species} is not supported")
data[:, simion_params.index("X")] = P.z * 1e3
data[:, simion_params.index("Y")] = P.y * 1e3
data[:, simion_params.index("Z")] = -P.x * 1e3
px = P.pz
py = P.py
pz = -P.px
data[:, simion_params.index("KE")] = P.kinetic_energy # [eV]
data[:, simion_params.index("AZ")] = np.arctan2(-pz, px) * (180 / np.pi) # [deg]
data[:, simion_params.index("EL")] = np.arctan2(py, np.sqrt(px**2 + pz**2)) * (
180 / np.pi
) # [deg]
data[:, simion_params.index("CWF")] = (
P.weight / e_
) # Charge Weighting Factor, derive from particle group weights
data[:, simion_params.index("COLOR")] = np.full(len(P), color)
# fname, X, fmt='%.18e', delimiter=' '
np.savetxt(filename, data, delimiter=",", header=header, comments="", fmt=" %.9e")
from scipy.constants import physical_constants
mc2 = 1e6 * physical_constants["electron mass energy equivalent in MeV"][0]
e_ = physical_constants["elementary charge"][0]
me = physical_constants["electron mass in u"][0]
def particle_group_to_SIMION(P, filename, color=0):
header = ";0"
simion_params = [
"TOB",
"MASS",
"CHARGE",
"X",
"Y",
"Z",
"AZ",
"EL",
"KE",
"CWF",
"COLOR",
]
# simion_units = {
# "TOB": "usec",
# "MASS": "amu",
# "CHARGE": "e",
# "X": "mm",
# "Y": "mm",
# "Z": "mm",
# "AZ": "deg",
# "EL": "deg",
# "CWF": "",
# "COLOR": "",
# }
data = np.zeros((len(P), len(simion_params)))
data[:, simion_params.index("TOB")] = P.t * 1e6 # [P.t] = sec, convert to usec
if P.species == "electron":
data[:, simion_params.index("MASS")] = np.full(len(P), me)
data[:, simion_params.index("CHARGE")] = np.full(len(P), -1)
else:
raise ValueError(f"Species {P.species} is not supported")
data[:, simion_params.index("X")] = P.z * 1e3
data[:, simion_params.index("Y")] = P.y * 1e3
data[:, simion_params.index("Z")] = -P.x * 1e3
px = P.pz
py = P.py
pz = -P.px
data[:, simion_params.index("KE")] = P.kinetic_energy # [eV]
data[:, simion_params.index("AZ")] = np.arctan2(-pz, px) * (180 / np.pi) # [deg]
data[:, simion_params.index("EL")] = np.arctan2(py, np.sqrt(px**2 + pz**2)) * (
180 / np.pi
) # [deg]
data[:, simion_params.index("CWF")] = (
P.weight / e_
) # Charge Weighting Factor, derive from particle group weights
data[:, simion_params.index("COLOR")] = np.full(len(P), color)
# fname, X, fmt='%.18e', delimiter=' '
np.savetxt(filename, data, delimiter=",", header=header, comments="", fmt=" %.9e")
In [52]:
Copied!
particle_group_to_SIMION(P, "text_ion_file.ion")
particle_group_to_SIMION(P, "text_ion_file.ion")
In [53]:
Copied!
def read_simion_ION_file(filename):
data = np.loadtxt(filename, comments=";", delimiter=",", skiprows=1)
simion_params = [
"TOB",
"MASS",
"CHARGE",
"X",
"Y",
"Z",
"AZ",
"EL",
"KE",
"CWF",
"COLOR",
]
return {simion_params[ii]: data[:, ii] for ii, p in enumerate(simion_params)}
def read_simion_ION_file(filename):
data = np.loadtxt(filename, comments=";", delimiter=",", skiprows=1)
simion_params = [
"TOB",
"MASS",
"CHARGE",
"X",
"Y",
"Z",
"AZ",
"EL",
"KE",
"CWF",
"COLOR",
]
return {simion_params[ii]: data[:, ii] for ii, p in enumerate(simion_params)}
In [54]:
Copied!
ions1 = read_simion_ION_file("text_ion_file.ion")
ions2 = read_simion_ION_file("test_ion_writer.ion")
ions1 = read_simion_ION_file("text_ion_file.ion")
ions2 = read_simion_ION_file("test_ion_writer.ion")
In [55]:
Copied!
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
px, py, pz = P["px"], P["py"], P["pz"]
p = np.sqrt(px**2 + py**2 + pz**2)
theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2
x = px / p
y = py / p
z = pz / p
ax.scatter(x[::], y[::], z[::], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
px, py, pz = P["px"], P["py"], P["pz"]
p = np.sqrt(px**2 + py**2 + pz**2)
theta = np.arctan2(py, px)
phi = np.arccos(pz / p)
hist, pedges = np.histogram(phi, bins=100, density=True)
pcs = (pedges[1:] + pedges[:-1]) / 2
x = px / p
y = py / p
z = pz / p
ax.scatter(x[::], y[::], z[::], ".")
ax.set_xlabel(r"$\hat{p}_x$")
ax.set_ylabel(r"$\hat{p}_y$")
ax.set_zlabel(r"$\hat{p}_z$")
Out[55]:
Text(0.5, 0, '$\\hat{p}_z$')
In [56]:
Copied!
for p in ions1:
print(f"{p}:", max(np.abs(ions1[p] - ions2[p])))
for p in ions1:
print(f"{p}:", max(np.abs(ions1[p] - ions2[p])))
TOB: 0.0 MASS: 0.0 CHARGE: 0.0 X: 0.0 Y: 0.0 Z: 0.0 AZ: 0.0 EL: 0.0 KE: 0.0 CWF: 0.0 COLOR: 0.0
In [57]:
Copied!
plt.plot(ions1["KE"], ions1["KE"] / ions2["KE"], ".")
plt.xlabel("KE (eV)");
plt.plot(ions1["KE"], ions1["KE"] / ions2["KE"], ".")
plt.xlabel("KE (eV)");
In [58]:
Copied!
os.remove("text_ion_file.ion")
os.remove("KEdist.txt")
os.remove("test_ion_writer.ion")
os.remove("text_ion_file.ion")
os.remove("KEdist.txt")
os.remove("test_ion_writer.ion")
In [ ]:
Copied!
In [ ]:
Copied!